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ABSTRACT 

Large  space  structures  are  difficult  to  control  because  of  the  high  order  of  their 
mathematical  models.  The  high  order  mathematical  model  makes  the  use  of  a  reduced 
order  model  to  control  the  structure  desirable.  The  Karhunen-Loeve  expansion  along 
with  Galerkin's  method  is  used  to  generate  a  reduced  order  model.  A  control  algorithm 
is  achieved  by  applying  linear  quadratic  regulator  theory  to  the  reduced  order  model. 

The  Karhunen-Loeve  basis  functions  or  mode  shapes  must  first  be  found  to  identify 
the  reduced  order  model.  Previous  results  have  shown  that  in  the  limit  as  the  structural 
damping  approaches  zero  the  Karhunen-Loeve  mode  shapes  and  natural  mode  shapes 
converge.  Numerical  techniques  are  applied  to  evaluate  the  structural  damping  required 
for  convergence.  Once  the  Karhunen-Loeve  mode  shapes  are  determined,  the  reduced 
order  control  model  is  applied  to  the  full  order  system.  The  performance  of  various 
Karhunen-Loeve  models  is  compared  by  measuring  the  modal  energies  in  the  controlled 
and  uncontrolled  modes. 


Ul 


C4 

THESIS  DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may  not 
have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made,  within 
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I.     INTRODUCTION 

A.  LARGE  SPACE  STRUCTURES 

The  lightweight,  flexible  materials  used  to  construct  large  space  structures  (LSS), 
like  a  space  station,  are  lightly  damped  and  when  disturbed  will  vibrate  for  a  consider- 
able amount  of  time.  This  prolonged  vibration  could  jeopardize  the  structural  integrity 
of  the  structure  or  disturb  experiments  on  the  LSS.  The  purpose  of  this  thesis  is  to  study 
the  effects  of  controlling  this  structural  vibration  with  a  reduced  order  model,  specifically 
a  Karhunen-Loeve  reduced  order  model. 

B.  PROBLEM  APPROACH 

The  solution  of  the  vibration  control  problem  requires  the  use  of  a  mathematical 
model  that  describes  the  behavior  of  the  system  in  time.  The  space  structure  is  modeled 
as  a  combination  of  small  plates  of  unit  mass  connected  to  form  the  complete  structure. 
The  vibrational  motion  of  these  plates  can  be  modeled  as  a  set  of  coupled  damped  har- 
monic oscillators.  Using  modal  analysis  the  mathematical  model  of  the  structure  can 
be  decoupled  to  yield  a  set  of  uncoupled  simultaneous  second  order  differential 
equations. 

For  LSS  this  model  is  of  very'  high  order,  making  control  design  difficult.  It  is 
therefore  desirous  to  use  a  reduced  order  model  (ROM)  to  control  the  structure.  A 
control  system  is  designed  for  the  LSS  using  the  Karhunen-Loeve  (reduced  order)  model. 
The  control  system  is  then  applied  to  the  LSS  after  it  has  been  disturbed  by  an  impulse. 
The  performance  of  this  system  and  a  control  system  based  on  a  modal  model,  truncated 
by  frequency,  are  compared. 

The  LSS  used  in  this  thesis  is  a  dual  keel  space  station,  see  Figure  1  on  page  2.  A 
mathematical  model  of  this  space  station  is  provided  courtesy  of  McDonnel  Douglas 
Astronautics.  A  computer  simulation  of  this  space  station  is  used  to  determine  the  ef- 
fectiveness of  the  control  system.  Karhunen-Loeve  (reduced  order)  models  of  increasing 
size  are  simulated  and  the  modal  energies  are  calculated  and  used  to  determine  the  rela- 
tive effectiveness  of  the  control  models.  These  results  will  then  be  compared  to  results 
obtained  for  the  modal  model  truncated  by  frequency. 


Figure  1. 


Representation  of  a  Dual  Keel  Space  Station  from  [Ref.  1:  p.  1] 


C.     ORGANIZATION 

•  Chapter  II  is  a  statement  of  the  mathematical  model  of  the  space  station. 

•  Chapter  III  discusses  the  Karhunen— Loeve  expansion,  the  relationship  between 
the  Karhunen— Loeve  mode  shapes  and  the  natural  mode  shapes  and  how  a  re- 
duced order  model  can  be  synthesized  using  the  Karhunen— Loeve  expansion  and 
Galerkin's  method. 

•  Chapter  IV  discusses  the  determination  of  the  Karhunen-Loeve  modes. 

•  Chapter  V  describes  the  control  solution  and  presents  the  control  simulation  re- 
sults. 

•  Chapter  VI  discusses  the  conclusions  based  on  the  simulation  results  and  a  re- 
commended design  procedure. 


II.     THE  MATHEMATICAL  MODEL 

A.  INTRODUCTION 

Large  space  structures  are  flexible,  lightly  damped  structures  that  can  vibrate  for  a 
considerable  amount  of  time  when  disturbed  by  an  external  force.  To  effectively  control 
such  a  system,  a  mathematical  model  describing  the  evolution  of  that  system  in  time  is 
required. 

B.  THE  MODAL  MODEL 

The  space  station  is  modeled  as  a  combination  of  small  plates  of  unit  mass  con- 
nected to  form  the  complete  structure.  This  model  can  be  visualized  as  a  system  of 
masses  connected  by  springs  and  dashpots  (the  springs  representing  the  stiffness  factor 
and  the  dashpots  representing  the  damping  factor).  The  displacement  of  the  masses  can 
be  described  by  a  second  order  matrix  differential  equation  of  motion: 


d_ 


Mx\t)  +  -g-  Kx{t)  +  Kx{t)  =  F(t)  (1) 


where 

x  is  the  generalized  coordinate  vector 
M  is  the  diagonal  system  mass  matrix 

-rr-  K  is  the  structural  damping  term 

d  is  the  damping  constant 

M,  is  the  frequency  of  oscillation  of  the  system 

K  is  the  symmetric  system  stiffness  matrix 

F(t)  is  the  system  forcing  function 

This  equation  represents  a  system  of  simultaneous,  second  order  differential  equations 
that  are  coupled  by  the  K  matrix. [Ref.  1:  p.  3] 

This  equation  can  be  uncoupled  and  the  system  represented  by  a  set  of  independent 
second  order  differential  equations.  This  is  done  through  the  process  of  modal  analysis 
which  is  outlined  in  [Refs.  2,  3].  The  resulting  modal  model  consists  of  a  set  of  inde- 
pendent second  order  differential  equations: 

[ii  +  ciQ.i}  +  Q2>1  =  XTF~]  (2) 


where 

•  r\  is  the  coordinate  vector  or  modal  amplitude  vector 

•  Q2  =  diag[wol,coo2.  ...  ,  wj 

•  XT  =  [*]  x2    ...  xn~]  the  transpose  of  the  modal  matrix  or  mode  shape  vector 

•  F  is  the  torquing  force  applied  at  a  point 

Next  a  discrete-time  equation  describing  the  motion  of  the  space  station  in  terms  of  its 
natural  modes  of  vibration  is  developed.  The  discrete-time  state  equation  for  the  ith 
equation  of  motion  is: 

Z^kT  +  1)  =  OX7)  Z/Lkl)  +  W)  xf[F(kT)  +  W{kT>]  (3) 

where 

Z,  is  a  vector  of  the  ith  modal  amplitude  and 
the  ith  modal  velocity 

O,  is  the  ith  state  transition  matrix 

r,  is  the  ith  input  vector 

X,r  is  the  transpose  of  the  ith  mode  shape  vector 

F  is  the  control  torque  force  vector  applied  at  a  point 

T  is  the  sampling  time 

k  is  the  time  index 

W  is  the  disturbance  input 

This  equation  is  used  for  computer  simulation  of  the  space  station  and  control 
solution. [Ref.  1:  p.  4] 


III.     APPLICATION  OF  THE  KARHUNEN-LOEVE  EXPANSION  TO  THE 
REDUCED  ORDER  CONTROL  OF  LARGE  SPACE  STRUCTURES 

A.  INTRODUCTION 

Large  flexible  structures,  such  as  a  space  station,  as  a  class  of  distributed  parameter 
systems  (DPS)  require  a  finite  dimensional  model  for  control  design.  This  model  may 
be  achieved  by  approximating  the  state  of  the  LSS  using  the  Karhunen-Loeve  (KL)  ex- 
pansion. The  expansion  is  truncated  to  provide  the  finite  dimensional  approximation 
of  the  state  for  control  desisn.  The  KL  model  that  results  describes  the  evolution  of  the 
approximated  state  of  the  structure. 

The  natural  mode  shapes  are  normally  used  for  modeling  and  control  of  flexible 
structures.  The  relationship  between  the  natural  mode  shapes  and  the  KL  mode  shapes 
is  described  in  this  chapter  as  well  as  the  use  of  Galerkin's  method.  Galerkin's  method 
is  used  to  generate  a  reduced  order  model  (ROM),  using  both  the  natural  mode  shapes 
and  the  KL  expansion. 

B.  THE  KARHUNEN-LOEVE  EXPANSION 

The  purpose  of  the  KL  expansion  is  stated  by  Stark  and  Woods  [Ref.  4:  p.  322]. 
"The  idea  [of  the  KL  expansion]  is  to  decompose  a  general  second-order  random  proc- 
ess into  an  orthonormal  expansion  whose  coefficients  are  uncorrected  random 
variables."  The  state  of  a  large  space  structure  (LSS),  y(x),  can  be  modeled  as  random 
process  since  it  depends  on  random  excitations,  i.e.,  noise  from  onboard  machinery  and 
actuators.  The  second  order  moments  of  the  LSS  are  proportional  to  the  physical  en- 
ergy and  therefore  exist.  The  LSS  can  therefore  be  approximated  using  the  KL  expan- 
sion [Ref.  5:  p.  12].  This  is  done  by  projecting  the  random  process  onto  an  orthonormal 
basis  and  truncating  to  N  terms.  The  value  chosen  for  N  is  a  matter  of  "engineering 
judgment".   [Ref.  5:  p.  13] 

The  selection  of  the  orthonormal  basis  is  made  by  solving  the  KL  eigenequation: 

<  Ryy{x,z),  4>fa)  >  =  k&fac)  (4) 

where 

•  4>,(x)  is  referred  to  as  an  eigenfunction  (or  the  KL  mode  shapes) 

•  X,  is  the  eigenvalue  and  is  a  measure  of  the  excitation  of  the  ith  basis  function 


•  Ryy(x,z)  is  the  correlation  function  of  y(x):  Ry})(x,z)  =  E[y{x)  y^x)'] 

•  <  • ,  •  >  is  an  inner  product:    <  a(z),  b(z)  >  =  f  aT(z)  M(z)  b(z)  dz 

•  M(z)  is  the  mass  density  of  the  structure 

•  Q.  is  the  spatial  extent  of  the  structure 


The  state  of  the  LSS  is  approximated  by: 


N 

ya  -  Yj>Mx)  (5) 

l=\ 

where 

•  y„  is  an  approximation  to  the  state  of  the  space  structure 

•  £,  is  a  set  of  coordinates  found  by  c,  =  <  4>(x), y(x)  > 

•  (f),(x)  is  the  ith  basis  function  or  KL  mode  shapes 

The  expansion  is  truncated  by  keeping  the  eigenfunctions  (KL  mode  shapes)  associated 
with  the  N  largest  eigenvalues.   [Ref.  5:  p.  13] 

The  KL  expansion  yields  the  best  approximation  to  the  random  process,  i.e.,  mini- 
mizes the  expected  value  of  the  norm  of  the  error,  of  any  orthogonal  expansion.  The 
approximation  error  is  defined  as: 

oo 

e(x)  =  y(x)  -  ya(x)  =    £  Cj4M  (6) 

For  a  proof  of  the  optimality  of  the  KL  expansion  see  [Ref.  6:  p.  11]  or  [Ref.  5:  p.  15]. 

C.     RELATIONSHIP  BETWEEN  THE  KARHUNEN-LOEVE  MODE  SHAPES  AND 
THE  NATURAL  MODE  SHAPES 

The  relationship  between  the  KL  mode  shapes  or  KL  basis  functions  and  the  natural 
mode  shapes  is  taken  from  Burl  [Ref.  5:  p.  13].  The  KL  mode  shapes  of  a  structure  are 
related  to  the  natural  mode  shapes  of  that  structure  by  a  linear  transformation  repres- 
enting a  change  of  basis  which  can  be  written: 


<f>i{x)  = 


oo 

-MV 


+  c? 


0 

Vi{x) 


=  VT(x)ci  (7) 


where 


v\*)  = 


n\{x)      G      niix)      G      i3{x) 

0        titix)        0        y2{.x)        0 


T      r   1        2        3 


r     r 


] 


(8) 


(9) 


and  {r},{x)}  is  the  set  of  natural  mode  shapes.  The  state  of  a  structure  consists  of  a 
generalized  position  and  velocity  which  can  be  expanded  in  terms  of  the  natural  mode 
shapes: 


y{x,t)  =  2jxp) 


j=\ 


n0 

0 


+  «/(0 


0 


=  »/r(x)a(/) 


(10) 


where 


r, 


a(/)  =  [a,(0     aj(0     a2(/)     a2(/)     a3(0      •••  ] 


(11) 


are  the  coordinates  and  velocities  of  the  natural  mode  shapes.    The  vectors,  c„  can  be 
found  by  solving  the  equation: 


£[KW/(r)]QC/  =  ;,C/ 


(12) 


where 


Q  =  diag 


"w?    0~ 

> 

2    rt 
cog     0 

_G      1 

j  *■•    » 

^56      G 

(13) 


Equations  10  and  11  give  the  KL  mode  shapes  in  terms  of  the  natural  mode  shapes. 
Equation  12  is  an  infinite  dimensional  eigenvalue  problem  that  can  be  solved  practically 
by  truncating  it  to  the  most  significant  terms. 

D.     REDUCED  ORDER  MODELING  USING  THE  KARHUNEN-LOEVE  MODES 

Galerkin's  method  is  used  to  produce  reduced  order  state  equations.    This  can  be 
done  with  either  the  KL  modes  or  the  natural  modes. 


The  discrete-time  state  equation  of  a  distributed  parameter  system  is: 

y(x,k)  =  Fy(x,k-\)  +  GJ{k)  (14) 

where 

y(x,k)  is  the  state 

f(k)  is  the  input 

k  is  the  time  index 

F  is  an  operator  on  the  state  space 

G  is  an  operator  from  0tm  to  the  state  space 

For  each  k,  y(x,k)  e  the  state  space,  f(k)  e  Skm 

A  finite  dimensional  approximation  to  this  equation  can  be  obtained  using  Galerkin's 
method 

Pny(x,k)  =  (PnFPn)y(x,k  -\)  +  PnG  u(k)  (15) 

where  P„  can  be  written  in  terms  of  a  basis  \_r\,  \  i  =  1,2, ...  ,  n~] 

n 

Pn(-)  =  YJVi(x)<  nfr),   •>  (16) 

The  KL  mode  shapes  (basis  functions)  can  be  used  in  the  above  equations  and  the  KL 
model  results.  The  natural  mode  shapes  can  be  used  in  the  above  equations,  which  is 
equivalent  to  truncating  the  modal  equations,  producing  the  modal  model. [Ref.  5  :  p. 
13] 

E.    SUMMARY 

The  KL  expansion  can  be  used  to  approximate  the  state  of  a  LSS.  This  is  done  by 
projection  of  the  state,  which  is  modeled  as  a  random  process,  onto  an  orthonormal 
basis  function.  The  basis  function  can  be  found  by  solving  the  KL  eigenequation.  These 
KL  basis  functions  will  yield  the  best  approximation  to  the  state  of  the  LSS.  Then  using 
Galerkin's  method  a  reduced  order  model  of  the  LSS  is  produced.  This  reduced  order 
model  is  used  to  generate  a  control  which  is  applied  to  the  entire  system. 


IV.     DETERMINATION  OF  THE  KARHUNEN-LOEVE  BASIS 

FUNCTIONS 

A.  INTRODUCTION 

The  KL  basis  functions  (or  mode  shapes)  can  be  determined  by  solving  the 
eigenequation  12.  This  is  a  tedious  and  laborious  processes.  It  can  be  shown  that  for 
lightly  damped  structures  the  KL  mode  shapes  can  be  determined  from  the  open  loop 
response  by  ordering  the  natural  mode  shapes  in  order  of  decreasing  modal  energies. 

The  KL  mode  shapes  were  calculated  using  the  KL  mode  program  in  Appendix  B. 
These  calculations  were  compared  to  the  mode  shapes  selected  by  observing  the  open 
loop  response  to  verify  that  this  is  a  valid  method  of  determining  the  KL  mode  shapes. 

The  value  of  the  structure's  damping  factor  used  in  the  program  was  increased  until 
the  KL  mode  shapes  no  longer  were  the  same  as  those  determined  from  the  open  loop 
response.  This  was  done  to  determine  how  large  the  structure's  damping  factor  could 
be  and  still  have  the  KL  mode  shapes  converge  to  the  natural  mode  shapes. 

B.  NUMERICAL  DETERMINATION  OF  THE  KARHUNEN-LOEVE  MODE 
SHAPES 

The  KL  mode  shapes  are  found  by  determining  the  eigenvalues  and  eigenvectors 
(equation  12)  of  the  covariance  matrix  £[a,(0  a/OX  This  is  an  infinite  dimensional 
matrix  which  is  truncated  to  a  finite  dimensional  square  matrix.  There  are  three  terms 
in  E[_o.{t)  ot.T{t)~\,  where  «(/)  is  defined  by  equation  11.  The  first  of  these  terms  is  com- 
puted, for  a  white  noise  input 


EC 


Too 
«Kt)«/t)]=  hKT)h/T)3T  (17) 


where  h,(t)  is  the  impulse  response  for  a,.(t)  [Ref.  6:  p.  66].  The  impulse  response  is  given 
by; 

Ht)  = %=j     e-^''   sinK/x/l-y2t)  (18) 


where 


•  b2,  is  the  modal  slope  in  the  x  direction 

•  y  is  the  damping  coefficient;  y  =  —  is  a  constant 

•  o)0,  is  the  natural  frequency 

Performing  the  integration  in  equation  17  produces  the  result: 


E[a/(t)  cc/t)]  =  Kf  K, 


a 


L  2[a2  +  b2]        2[a2  +  c2]   J 


where 


•    K,„  = 


b2l,j 


•    a  =  y  ojoi  +  y  cocj 


•  b  =  (D0lJ\-y2  -cojl-y2 

•  c  =  cociyJ\  -y2  +toorJl  -y2 

The  other  terms  in  ECc<(t)  otr(t)~J  can  be  found  in  a  similar  way.   They  are: 


ECa,.(t)  oc/t)]  =  K,  K;  Li 


+ 


L  2[a2  +  b2;         2[a2  +  c2] 
-  K;  Kj  M, 


L  2[a2  +  b2]        2[a2  +  c2]   J 


and 


E[dl<t)«Xt)]  =  K/K/ 


L,L,a  -  L,M  b  -  LjM,b  +  M.M.a 


+ 


'i^r 


2[a2  +  b2] 


L,L,a-L;M,c-M;M;a 


i-^r 


2[a2  +  c2] 


where 


(19) 


(20) 


(21) 


•  LtJ  =  Q)eUrJl-y2 

•  MiJ  =  ya>eitJ 

Equations  19,  20,  and  21  specify  the  eigenvalue  problem  equation  12. 

The  eigensolution  specifies  the  transformation  from  the  natural  mode  shapes  to  the 
KL  mode  shapes.  This  finite  dimensional  eigenvalue  problem  is  solved  and  the  KL  mode 
shapes  are  approximated  as  a  linear  combination  of  the  first  fifty  (flexible)  natural  mode 
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shapes.     The  program  in  Appendix  B  computes  the  KL  mode  shapes  by  solving  the 
eigenvalue  problem  equation  12. 

C.     EMPIRICAL  DETERMINATION  OF  THE  KARHUNEN-LOEVE  MODE 
SHAPES 

The  numerical  computation  of  the  KL  mode  shapes  is  difficult.  Solving  the 
eigenvalue  problem  requires  the  determination  of  the  impulse  response  and  solving  for 
all  the  terms  in  the  covariance  matrix  E[a,(t)  a,(t)].  It  turns  out  that,  in  the  limit  as  the 
damping  factor  of  the  structure  approaches  zero,  the  KL  mode  shapes  converge  to  the 
natural  mode  shapes  [Ref.  7j. 

If  the  modal  energies  for  the  open  loop  response  are  computed  and  the  mode  shapes 
ordered  by  decreasing  modal  energy,  these  mode  shapes  should  be  the  same  as  those 
computed  numerically.   The  output  of  the  system,  y,  is  defined: 


yi 


>56 


W7      0 

0       1 


0      0 
0      0 


CO 


0 

0 

0 

0 

'*7~ 

• 
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• 

• 

• 

• 

• 

• 

• 

J56 

0 

■*56 

0 

1 

(22) 


where  yt  e  0l2  and 


£[ith  Modal  Energy]  =  £[  ||j/,  ||2]. 


(23) 


The  energy,  given  in  equation  23  can  be  written  [Ref.    8  ]: 


£[IU 


1  II2]  -   f  Ify 


(t)fdt 


(24) 


where  hy,  is  the  response  of  y,  due  to  an  impulse  applied  at  the  disturbance  input  (node 
69  or  55).  Equation  24  is  evaluated  using  computer  simulation.  The  mode  shapes  found 
are  compared  to  those  determined  numerically  for  noise  input  at  node  69  or  node  55,  see 
Table  1  on  page  12. 
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Table   1.       DETERMINATION       OF       THE 
KARHUNEN-LOEVE  MODE  SHAPES 


FIRST       FIFTY       FLEXIBLE 


Karhunen-Loeve 
mode  shapes  deter- 
mined from  the  open 
loop  response. 

Karhunen-Loeve 
mode  shapes  numer- 
ically calculated. 

Noise  input  Location: 

Node  55 

Node  69 

Node  55 

Node  69 

40 

54 

40 

54 

43 

51 

43 

51 

■? 

31 

7 

31 

17 

52 

17 

52 

28 

36 

28 

36 

15 

7 

15 

7 

44 

30 

44 

30 

35 

28 

35 

28 

39 

48 

39 

48 

33 

26 

33 

26 

41 

35 

41 

35 

31 

50 

31 

50 

42 

15 

42 

15 

45 

55 

45 

55 

53 

41 

53 

41 

48 

56 

48 

56 

11 

44 

11 

44 

36 

53 

36 

53 

26 

33 

26 

}}, 

21 

34 

21 

34 

8 

25 

8 

25 

16 

37 

16 

37 

38 

38 

38 

38 

29 

46 

29 

46 

50 

27 

50 

27 

30 

43 

30 

43 

23 

40 

23 

40 

19 

23 

19 

23 

32 

21 

32 

21 

25 

16 

25 

16 

51 

8 

51 

8 

27 

45 

27 

45 

22 

42 

25 

42 

37 

11 

37 

11 

20 

47 

20 

47 

10 

49 

10 

49 

52 

18 

52 

18 

55 

20 

55 

20 

9 

29 

46 

29 

46 

10 

49 

10 

49 

24 

9 

24 

34 

17 

34 

17 

56 

39 

56 

39 

24 

22 

24 

22 

54 

32 

54 

32 

13 

9 

13 

9 

18 

13 

18 

13 

47 

14 

47 

14 

14 

19 

14 

19 

12 

12 

12 

12 

12 


Table  1  shows  that  the  mode  shapes  are  the  same.  Therefore,  the  KL  mode  shapes  can 
be  determined  from  the  open  loop  response  when  the  structural  damping  factor  is  suffi- 
ciently small.  In  the  case  of  the  space  station  simulated  in  this  thesis  a  value  of  0.001 
was  used.  The  next  question  is.  how  big  can  the  damping  factor  of  the  structure  be  be- 
fore the  KL  mode  shapes  fail  to  converge  to  the  natural  mode  shapes? 
1.     Required  Magnitude  of  the  Damping  Factor  for  Convergence 

When  the  structural  damping  factor  is  small  enough,  the  eigenvectors  solved  for 
in  equation  12  approximate  the  natural  basis 

e|  =  n  1  0  0  ...  ]  (25) 

or 

eJ  =  rjO  1  0  0  ...  ]  (26) 

etc.;  one  element  is  unity  and  the  other  elements  are  zero.  When  these  vectors  are  sub- 
stituted into  equation  7  it  is  easy  to  see  that  the  KL  mode  shapes  are  the  natural  mode 
shapes.  If  the  damping  factor  is  increased,  at  some  value  the  eigenvectors  have  more 
than  one  non-zero  element  and  the  KL  mode  shapes  become  a  linear  combination  of 
natural  mode  shapes. 

The  damping  factor  was  increased  successively  by  a  factor  of  two  from  a  starting 
value  of  0.0005.  The  norm  of  the  error  (or  error  norm)  is  used  as  a  means  of  measuring 
how  closely  the  eigenvectors  obtained  from  the  KL  mode  program  approach  the  ideal 
of  equation  25.   The  error  norm  is  defined  as  follows: 


50 

k  =  minimumll  ck  —  ej  ||  =     /(l  —  Xj)2  4-  y  x2  (27) 

Mi 

/=2 


where 


ck  =  [xi  *2  x3  -  xj  -  xso]  (28) 

and  for  simplicity  x;  is  assumed  to  be  positive  and  the  largest  component  of  c,.  Equation 
27  can  be  simplified  to 
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£||k=    /l--2x.+    >  x?  (29) 


Since  the  eieenvectors  are  normalized, 


50 

y x?  =  i  (30) 


and  substituting  equation  30  into  equation  29  yields 


||£||k  =  x/2(l  -xp  (31) 

a.  Results 

The  norm  of  the  error,  £,  and  the  natural  mode  shape  sequence  obtained 
from  the  KL  mode  program  were  used  to  determine  how  large  the  damping  factor  could 
be  and  still  have  the  KL  mode  shapes  and  natural  mode  shapes  converge.  The  error 
norm  is  presented  graphically  for  each  KL  mode;  all  cases  are  for  disturbance  torques 
due  to  actuator  noise.   The  cases  presented  are: 

•  Figure  2  on  page  16;  for  a  damping  factor  of  0.0005  all  but  two  KL  modes  have 
an  error  norm  below  0.2.  The  natural  mode  sequence  is  as  in  Table  1  on  page  12. 

•  Figure  3  on  page  17;  for  a  damping  factor  of  0.001  nearly  74  percent  of  the  KL 
modes  have  an  error  norm  of  0.2  or  below,  the  natural  mode  sequence  is  presented 
in  Table  1  on  page  12.  The  natural  mode  shapes  are  a  good  approximation  of  the 
KL  mode  shapes. 

•  Figure  4  on  page  18;  for  a  damping  factor  of  0.002  only  46  percent  of  the  KL 
modes  have  an  error  norm  at  or  below  0.2,  and  the  natural  mode  sequence  no 
longer  conforms  to  that  in  Table  1  on  page  12,  i.e.,  the  KL  modes  and  natural 
modes  are  not  conversing. 

•  Figure  5  on  page  19;  for  a  damping  factor  of  0.005  only  28  percent  of  the  KL 
modes  have  an  error  norm  of  0.2  or  less  and  the  KL  modes  and  the  natural  modes 
are  more  divergent. 

b.  Conclusion 

If  the  damping  factor  is  greater  than  0.001  the  natural  mode  shapes  are  not 
a  good  approximation  of  the  KL  mode  shapes.  The  convergence  criteria  is  that  the 
norm  of  the  error,  for  50  percent  of  the  KL  modes  or  greater,  is  more  than  0.2.  The 
value  for  the  damping  factor  used  in  the  space  station  simulations  as  noted  in  previous 
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chapters  was  0.001.  The  assumption  that  the  KL  mode  shapes  and  natural  mode  shapes 
converge  is  valid  for  this  damping  factor  and  determining  the  KL  modes  from  the  modal 
energies  is  a  good  approximation. 
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ERROR  NORM  VS  DAMPING  FACTOR 
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Figure  2.       Norm  of  the  error  values  for  the  KL  modes,  excitation  at  node  69  , 
damping  factor  =  0.0005 
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ERROR  NORM  VS  DAMPING  FACTOR 
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Figure  3.      Norm  of  the  error  values  for  the  KL  modes,  excitation  at  node  69  , 
damping  factor  =  0.001 
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ERROR  NORM  VS  DAMPING  FACTOR 
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Figure  4.       Norm  of  the  error  values  for  the  KL  modes,  excitation  at  node  69  , 
damping  factor  =  0.002 
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ERROR  NORM  VS  DAMPING  FACTOR 

o 

oi 

CO 

LEGEND 

1.4         1.6 

I              I 

DAMPING=0.005 

CN 

2    -" 

cr 

o 

Z  9_ 

cr  ' 

o 

q:  go 

CC  6- 

LU 

6~  _ 

1      r1    nn 

i    w 

T+- 

i     I 

i  n 

d~ 

pi 

Lrun 

Lr 

j 

(N 

-j 

d  - 

'— ' 

i— i 

Ln 

o 

_ 

— 

d 

51 

i       i       i      l       l      l      l      l      I      I       l      l      l       l       I       i 
0     3    6     9    12   15  18  2124  27  30  33  36  39  42  45  48 

MODE 

Figure  5.      Norm  of  the  error  values  for  the  KL  modes,  excitation  at  node  69  , 
damping  factor  =  0.005 
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V.     CONTROL  SOLUTION 

A.     INTRODUCTION 

The  response  of  the  space  station  to  disturbance  and  control  inputs  is  simulated 
using  the  modal  model  discussed  in  Chapter  II.  The  Fortran  program  that  simulates  this 
model  was  written  by  Preston  [Ref.  1:  p.  47],  and  is  used  in  this  thesis  with  minor  mod- 
ification, (see  Appendix  A).  The  objective  of  the  simulation  is  to  determine  the  system 
response  to  disturbances  applied  at  modes  69  and  55  (see  Figure  6)  on  the  structure, 
using  increasing  sizes  of  KL  models  (i.e.,  5,  10,  20  modes).  The  response  of  the  system 
is  depicted  graphically  by  displaying  the  energy  in  each  mode  in  english  units  of  inch- 
pounds  (in— lbs). 


Figure  6.      Disturbance  Locations 
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B.     FULL  ORDER  SPACE  STATION  MODEL 

The  plant  being  controlled  in  this  thesis,  as  stated  in  Chapter  I,  is  a  preliminary 
version  of  the  NASA  dual  keel  space  station.  NASTRAN  was  used  to  generate  the  first 
fifty  flexible  modes  (starting  with  mode  seven).  The  modal  model  that  forms  the  full 
order  model  is  composed  of  these  flexible  modes: 
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(32) 


where 


*/  = 


modal  amplitude 
modal  velocitv 


(33) 


and 


co,  is  the  natural  frequency  of  the  ith  mode 

d  =  0.001  is  the  damping  coefficient 

uef  is  the  control  input  torques  from  three  orthogonally  oriented  control  mo- 
ment gyros 

w  e  ^?3  is  a  random,  white  noise  disturbance  input  (three  orthogonal  disturbance 
torques) 

<j),{p)  e  0P  is  the  modal  slopes  at  p 

pc  is  the  location  of  the  control  moment  gyros  (node  69  on  Figure  6  on  page  20) 

pn  is  the  location  of  the  disturbance  input  (either  the  control  moment  gyro  location, 
node  69,  or  the  alpha  joint,  node  55,  as  seen  in  Figure  6  on  page  20). 
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This  model  is  used  as  the  full  order  model  in  the  simulation  of  the  space  station.  How- 
ever, it  is  itself  a  reduced  order  model  since  the  "infinite"  number  of  vibrational  modes 
of  the  structure  are  truncated  to  fifty. 

It  is  assumed  that  perfect  modal  amplitude  and  modal  velocity  information  is  avail- 
able. This  assumption  of  perfect  sensor  information  is  not  realistic,  but  isolates  the  ef- 
fect of  modal  truncation  on  control  algorithm  synthesis. 

C.  THE  REDUCED  ORDER  MODEL 

The  reduced  order  model  of  the  space  station  is  obtained  by  truncating  modes.  The 
criteria  for  truncating  modes  is  either  the  modal  frequency  or  the  Karhunen-Loeve  or- 
dering (the  KL  ordering  found  in  the  previous  chapter  to  be  equivalent  to  an  ordering 
based  on  the  energy  measured  during  open  loop  excitation). 

D.  PERFORMANCE  FUNCTION  AND  OPTIMAL  CONTROL 

The  full-order  model  and  a  reduced  order  model  have  been  established.  Next,  a 
mathematical  expression  of  system  performance  is  required.  The  performance  function, 
J,  stated  here  is  developed  in  detail  in  [Ref.  1:  p.  10  ].  J  consists  of  the  total  energy  in 
the  modeled  modes  plus  a  control  energy  term: 


J = [*T  •  •  •  */L*] 
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(34) 


where 

•  imax  is  the  number  of  modes  in  the  reduced  order  model 

•  r  =  10-12  is  a  weighting  coefficient  on  the  control  energy  term  selected  to  yield  time 
constants  on  controlled  modes  of  approximately  30  seconds. 

Vibration  damping  is  achieved  by  application  of  steady  state,  linear  quadratic  regulator 
theory  [Ref.  9]  to  the  KL  (reduced  order)  model.  The  control  torque  vector,  u(k),  is  the 
product  of  an  optimal  gain  matrix,  L,  and  the  time  varying  state  matrix,  Z  (equation  3). 
The  L  matrix  is  found  by  solution  of  the  Ricatti  equations  to  minimize  the  performance 
function,  J,  [Ref.  1:  p.  12). 
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E.    SIMULATION  RESULTS 

The  space  station  was  simulated  and  the  modal  energies  calculated  for  disturbance 
inputs  applied  at  nodes  69  and  55  as  shown  in  Figure  6  on  page  20.  The  control  system 
consists  of  KL  (reduced  order)  models  of  five,  ten  or  twenty  modes.  The  data  is  pre- 
sented graphically  as  energy  in  each  of  the  modes.  The  open  loop  response  is  included 
for  comparison  except  where  the  difference  in  scale  precludes  it.  The  closed  loop  re- 
sponse for  a  reduced  order  model  generated  by  modal  truncation  is  included  for  com- 
parison with  the  KL  model.   The  cases  presented  are: 

•  Open  loop  (no  control)  response  which  identifies  the  reduced  order  model  to  be 
used.    See  Figure  7  and  Figure  8. 

•  Closed  loop  response  for  a  reduced  order  model  generated  by  modal  truncation 
with  the  modes  ordered  by  natural  frequency.    See  Figure  9  to  Figure  14. 

•  Closed  loop  response  for  a  Karhunen-Loeve  (reduced  order)  model.  The  KL  modes 
were  determined  from  the  open  loop  responses.  These  results  are  shown  in 
Fisure  15  through  Fisure  19. 

•  Closed  loop  response,  with  a  KL  model  based  on  the  open  loop  response,  for  the 
system  excited  by  the  control  actuators.  These  results  are  shown  in  Figure  20 
through  Figure  22. 

These  results  are  discussed  in  detail  in  Chapter  VI. 
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Figure  7.      Open  Loop  Response,  excitation  at  node  55. 
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Figure  8.      Open  Loop  Response,excitation  at  node  69 
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Figure  9.      Five  Controlled  modes  selected  by  frequency,  excitation  at  node  55. 
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Figure  10.      Ten  controlled  modes  selected  by  frequency,  excitation  at  node  55. 
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Fi°ure  11.      Twenty  controlled  modes  selected  by  frequency,excitation  at  node  55. 
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Figure  12.      Five  controlled  modes  selected  by  frequency,  excitation  at  node  69. 
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Figure  13.      Ten  controlled  modes  selected  by  frequency,  excitation  at  node  69. 
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Figure  14.      Twenty  controlled  modes  selected  by  frequency,  excitation  at  node  69. 
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Figure  15.      Five  controlled  modes  selected  by  modal  energy,  excitation  at  node  55. 
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Figure  16.      Ten  controlled  modes  selected  by  modal  energy,  excitation  at  node  55. 
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Figure  17.      Twenty  controlled  modes  selected  by  modal  energy,  excitation  at  node 

55. 
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Figure  18.      Five  and  ten  controlled  modes  selected  by  modal  energy,  excitation  at 
node  69. 
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Figure  19.      Ten  and  twenty  controlled  modes  selected  by  modal  energy,  excitation 
at  node  69. 
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Figure  20.      Five  controlled  modes  selected  by  modal  energy  due  to  excitation  at  node 
69,  actual  excitation  at  node  55. 
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Figure  21.      Ten  controlled  modes  selected  by  modal  energy  due  to  excitation  at  node 
69,  actual  excitation  at  node  55. 
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Figure  22.      Twenty  controlled  modes  selected  by  modal  energy  due  to  excitation  at 
node  69,  actual  excitation  at  node  55. 
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VI.     CONCLUSIONS  AND  RECOMMENDATIONS 

A.     CONCLUSIONS 

Results  supported  theoretical  claim  that  the  KL  mode  shapes  converge  to  the  na- 
tural mode  shapes.  For  the  space  station  example,  a  damping  factor  of  0.001  was  re- 
quired to  make  the  approximation  based  on  convergence  good. 

Preston  discovered,  see  [Ref.  1:  p.  46],  that  truncating  the  modes  based  on  natural 
frequency  yielded  poor  results.  Modes  with  large  coupling  to  the  control  system  are  not 
modeled  when  natural  frequency  is  the  method  of  truncation.  Instead  modes  with  very 
little  coupling  to  the  control  system  are  modeled  and  in  attempting  to  control  these 
modes  very  large  control  torques  are  generated  causing  large  excitations  in  the  strongly 
coupled  but  unmodeled  modes  as  seen  in  Figures  9  through  14. 

The  excitation  of  each  mode  depends  on  two  factors:  the  natural  frequency  which 
determines  the  damping  and  the  amplitude  of  the  mode  shape  which  determines  how 
much  excitation  is  received  by  each  mode.  The  amplitude  of  the  mode  shape  is  the 
dominate  factor  which  is  evident  from  the  open  loop  responses  see  Figures  7  and  8. 

Truncating  based  on  the  open  loop  response,  i.e.,  the  Karhunen-Loeve  model  also 
yields  poor  results  except  when  it  is  developed  from  the  open  loop  response  for  the  case 
of  the  disturbance  torques  being  due  to  actuator  noise,  see  Figure  8  on  page  25.  When 
truncating  modes  based  on  the  open  loop  response,  the  coupling  from  the  noise  input 
is  the  dominate  factor  in  selecting  the  modes  to  include  in  the  model.  In  this  case  modes 
with  large  open  loop  excitation  are  controlled,  but  the  control  coupling  still  dominates 
yielding  large  excitation  in  the  unmodeled  modes  (see  Figures  15  through  17).  When  the 
disturbance  torques  are  due  to  actuator  noise,  the  modes  with  the  largest  open  loop 
excitation  are  also  the  modes  with  the  largest  control  coupling  and  the  control  system 
works  well  (see  Figures  18  and  19). 

The  KL  model  based  on  the  open  loop  response  for  disturbance  torques  due  to 
actuator  noise  works  well  when  applied  to  the  space  station  with  a  disturbance  torque 
applied  at  another  location.  As  seen  in  Figures  20  through  22,  the  control  system  drives 
the  modeled  modes  close  to  zero  without  exciting  other  modes  to  a  significant  degree. 
The  problem  with  this  configuration  is  that  modes  with  large  open  loop  excitation  are 
not  necessarily  modeled  and  may  be  left  excited  by  the  control  svstem. 
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B.     DESIGN  PROCEDURE 

This  design  procedure  is  based  on  the  above  conclusions  and  is  recommended  for 
controlling  vibration  in  LSS.  First,  the  reduced  order  control  model  should  be  based  on 
the  KL  mode  shapes  determined  from  the  open  loop  response  for  disturbance  torques 
due  to  actuator  noise.  Next,  check  the  open  loop  response  for  disturbance  torques  at 
locations  on  the  structure  where  noise  inputs  are  most  likely,  e.g.,  the  alpha  joint. 
Modes  strongly  coupled  to  the  disturbance  input  must  be  included  in  the  KL  model. 
This  will  determine  the  size  of  the  KL  (reduced  order)  model.  If,  after  all  this  is  done, 
the  results  are  still  unacceptable,  then  adding  additional  control  elements  should  be 
considered.  The  placement  of  the  additional  control  element(s)  is  determined  by  the 
node(s)  with  the  largest  modal  amplitude(s)  for  the  modes  that  remain  excited  by  the 
initial  control  system.   The  process  is  repeated  until  an  acceptable  control  is  achieved. 
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APPENDIX  A.     SPACE  STRUCTURE  SIMULATION  PROGRAM 

This  program  simulates  the  dual  keel  space  station  described  in  Chapter  I  by  im- 
plementing the  model  described  in  Chapter  IV.  The  control  described  in  Chapter  IV  is 
also  simulated.   For  a  detailed  explanation  of  this  program  and  its  development  see  [Ref. 
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C      ************************************  ft*  A****************  ***"A 

C 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


LAMA  =  VECTOR  OF  THE  SQUARE  OF  THE  NATURAL  FREQUENCIES 

UGVEX  =  MODE  POSITONS  AND  SLOPES  OF  THE  NODAL  POINTS 

PHI  =  STATE  TRANSITION  MATRICIES  FOR  EACH  MODE 

GAMMA  =  INPUT  TRANSITION  MATRIX 

A  =  DIAGONAL  MATRIX  CONSISTING  OF  PHI 

B  =  INPUT  MATRIX  OF  GAMMA  AND  CONTROL  NODE  SLOPES 

BN  =  NOISE  INPUT  MATRIX  OF  GAMMA  AND  NOISE  NODE  SLOPES 

DAMP  =  DAMPING  FACTOR 

SAMPT  =  SAMPLING  TIME 

IMPLSE  =  IMPULSE  INPUT  FUNCTION 

TCX,  TCY,  TCZ  =  CONTROL  TORQUE  VALUES 

IMPX,  IMPY,  IMPZ  =  AXIS  IMPULSE  NOISE  VALUES 

ENERGY  =  SYSTEM  ENERGY  COST  VALUE  FOR  A  GIVEN  POINT  IN  TIME 

CNTCST  =  SYSTEM  CONTROL  COST  VALUE  FOR  A  GIVEN  POINT  IN  TIME 

COST  =  TOTAL  SYSTEM  COST  VALUE  FOR  A  GIVEN  POINT  IN  TIME 

TOTCST  =  SYSTEM  COST  SUMMED  OVER  ALL  TIME 

MIN  =  NUMBER  OF  MINUTES  SYSTEM  WILL  BE  OBSERVED 


********** 


SAMPLE  OF  SPACEN  EXEC  FILE 


************** 


THIS  FILE  MUST  BEGIN  IN  COLUMN  1  AND  RUN  WITH  THE  FOLLOWING 
SEQUENCE  FOR  THE  INITIAL  RUN  OF  THE  PROGRAM: 


FORTVS2  SPACEN 
SPACEN 


(COMPILES  PROGRAM) 
(LOADS  AND  RUNS  PROGRAM) 


SUBSEQUENT  PROGRAM  RUNS  CAN  ELIMINATE  "FORTVS2  SPACEN*'  IF  NO 
CHANGES  HAVE  BEEN  MADE  TO  THE  PROGRAM. 

CP  DEF  STOR  2M 

FI  4  DISK  KLAMA  OUTPUT  B  (PERM 

30  DISK  XI  OUTPUT  A  (RECFM  F  BLOCK  80  PERM 


FI 
FI 
FI 
FI 


31  DISK  MODENG  SPACEN  A  (RECFM  F 

32  DISK  TORQUE  OUTPUT  A  (RECFM  F 

33  DISK  ENERGY  OUTPUT  A  (RECFM  F 
FI  34  DISK  MDECST  OUTPUT  A  (RECFM  F 
FI  35  DISK  COUNT  INPUT  A  (RECFM  F 


BLOCK  80  PERM 
BLOCK  80  PERM 
BLOCK  80  PERM 
BLOCK  80  PERM 


BLOCK  80  PERM 
FI  40  DISK  UTILITY  DATA  A  (RECFM  F  BLOCK  80  PER 
FI  41  DISK  RUN  DATA  A  (RECFM  F  BLOCK  80  PERM 


FI  42  DISK  KUG69  OUTPUT  A 
FI  43  DISK  KUG23  OUTPUT  A 
FI  44  DISK  KUG55  OUTPUT  A 
LOAD  SPACEN 
START  *  NOXUFLOW 


(RECFM  F 
(RECFM  F 


BLOCK  80  PERM 
BLOCK  80  PERM 


(RECFM  F  BLOCK  80  PERM 


*********************************************************** 

PI  =  4.  0D0  *  ATAN(l.ODO) 

SAMPT  =  0.  0 

DAMP  =  0.  0 

MODAL  =  0 

IMPX  =  0.  0D0 

IMPY  =  0.  0D0 

IMPZ  =  0.  0D0 

NF  =  100 
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NG  =  100 
NH  =  100 
NZ  =  200 

Q       *********** 

C 

C     *****  NUMBER  OF  MINUTES  THE  SYSTEM  WILL  BE  OBSERVED  ***** 

C 

MIN  =  120.0 
C 

Q  *********************************************************** 

C     *****  SET  LENGTH  OF  MODAL  MODEL  ***** 

Q  *********************************************************** 

c 

MODAL  =56 
C 

Q       *********************************************************** 

C     *****  READ  LAMA  MATRIX  ***** 

Q       *********************************************************** 

c 

READ(4,1001)  NAM 
READ(4,1002)(LAMA(I),I=1,100) 
C 

Q       *********************************************************** 

C     *****  SCREEN  INTERACTION  ***** 

Q       *********************************************************** 

c 
c 

C     ************       STARTING  MODE  NUMBER      ************** 
C 

SMODE  =  7 
C 
C     ***************    NUMBER  OF  MODES  TO  SCAN    ************* 

C 

MODE  =5 

EMODE  =  SMODE  +  MODE  -  1 
C 

C     ************    NOISE  INPUT  POSITION      ***************** 
C 

NODE  =55 

AXIS  =  1 
C 
C     *************        r  MATRIX  VALUE       **************** 

C 

RM  =  1E-12 
C 

C     ********». 01  FOR  FULL  SAMPLING  TIME  °.05  FOR  REDUCED  ****** 
C 

SAMPT  =  0.  01 
C 
C     ***************        DAMPING  FACTOR        ************** 

DAMP  =  0.001D00 
C 

DO  75  I  =  1,100 

READ(42,1040)  (UG69( I ,K) ,K=1 ,3) 

READ(43,1040)  (UG23( I ,K) ,K=1 ,3) 
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READ(44,1040)  (UG55( I ,K) ,K=1 ,3) 
75    CONTINUE 
q  ***************        BEGIN  RUN 


*************** 


* 
C 


C 

c 
c 


DO  505  M  =  1,3 

WRITE  (41,700)  SMODE 

WRITE  (41,701)  MODE 

WRITE  (41,706)  EMODE 

WRITE  (41,702)  NODE 

WRITE  (41,703)  RM 

WRITE  (41,704)  SAMPT 

WRITE  (41,705)  DAMP 

WRITE  (41,707)  MIN 

WRITE  (41,708)  MODAL 

NOISE  AXIS  INPUT  AND  LOCATION 


******** 

IF(AXIS.EQ. 
IMPX  =  1 

ELSEIF(AXIS 
IMPY  =  1 

ELSEIF(AXIS 
IMPZ  =  1 

ELSEIF(AXIS 
IMPX  =  1 
IMPY  =  1 
IMPZ  =  1 

ENDIF 


********** 


1)THEN 
.  ODO/SAMPT 
.  EQ.  2) THEN 
. ODO/SAMPT 
. EQ. 3)THEN 
. ODO/SAMPT 
.  EQ.  4)  THEN 
. ODO/SAMPT 
. ODO/SAMPT 
. ODO/SAMPT 


C 

c 
c 


45 
40 
C 


65 
60 
C 


70 
C 
C 
C 


COUNT  =  0 
*********** 

DO 


INITIALIZE  MATRICIES 


**************** 


40  I  =  1,188 
DO  45  J  =  1,188 

PHI(I,J)  =0.0 

CONTINUE 
CONTINUE 

DO  60  I  =  1,188 
DO  65  J  =  1,3 

B(I,J)  =0.0 
BN(I,J)  =0.0 

CONTINUE 
CONTINUE 

DO  70  K  =  7,100 

X1(K)  =0.0 

X2(K)  =0.0 

MODEN(K)  =0.0 

RMODEN(K)  =0.0 
CONTINUE 

*********************************************************** 
*****  BEGIN  MAIN  PROGRAM  ***** 
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C     *****      ESTABLISH  A,  B  AND  B"NOISEM  MATRICIES      ***** 

Q       *********************************************************** 

c 

DO  600  I  =  SMODE, MODAL 

WN  =  DBLE(SQRT(LAMA(I))) 
GMA  =  DAMP*WN/2. 0 
EGT  =  DEXP(-GMA*SAMPT) 
Wl  =  DSQRT((WN**2)-(GMA**2)) 
C0SW1T  =  DC0S(W1*SAMPT) 
SINW1T  =  DSIN(W1*SAMPT) 
C 

IF(WN.  EQ.  0)THEN 

PHII(1,1,I)  =  EGT*C0SW1T 

PHII( 1,2,1)  =  SAMPT 

PHII(2,1,I)  =  0 

PHII(2,2,I)  =  EGT*C0SW1T 


ELSE 


GAMMA(1,I)  =  0 
GAMMA(2,I)  =  0 


PHII(1,1,I)  =  EGT*(C0SW1T  +  (GMA*(W1**(-1)))*SINW1T) 

PHII(1,2,I)  =  (W1**(-1))*EGT*SINW1T 

PHII(2,1,I)  =  -(WN**2)*(W1**(-1))*EGT*SINW1T 

PHII(2,2,I)  =  EGT*(C0SW1T  -  (GMA*(W1**( -1)))*SINW1T) 
C 

GAMMA(1,I)  =  (WN**(-2))*(1.0DO-EGT*(COSW1T+(GMA/W1)*SINW1T)) 

GAMMA(2,I)  =  (W1**(-1))*EGT*SINW1T 
C 

ENDIF 
C 

600   CONTINUE 
C 

V  =  1 
C 

DO  610  K  =  SMODE, MODAL 
C 

PHI(V,V)  =  PHII(1,1,K) 

PHI(V,V+1)  =  PHII(1,2,K) 

PHI(V+1,V)  =  PHII(2,1,K) 

PHI(V+1,V+1)  =  PHII(2,2,K) 
C 

B(V,1)  =  GAMMA(1,K)*DBLE(UG69(K,1)) 

B(V,2)  =  GAMMA(1,K)*DBLE(UG69(K,2)) 

B(V,3)  =  GAMMA(1,K)*DBLE(UG69(K,3)) 

B(V+1,1)  =  GAMMA(2,K)*DBLE(UG69(K,1)) 

B(V+1,2)  =  GAMMA(2,K)*DBLE(UG69(K,2)) 

B(V+1,3)  =  GAMMA(2,K)*DBLE(UG69(K,3)) 
C 

V  =  V+2 
C 

610   CONTINUE 
C 

DO  605  I  =  1,100 
UGVEX(I,1)  =0.0 
UGVEX(I,2)  =  0.  0 
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c 
c 


UGVEX(I,3)  =0.0 

IF(NODE.EQ.  23)  THEN 

UGVEX(I,1)  =  UG23(I,1) 
UGVEX(I,2)  =  UG23(I,2) 
UGVEX(I,3)  =  UG23(I,3) 
ELSEIF(NODE.EQ. 55)  THEN 
UGVEX(I,1)  =  UG55(I,1) 
UGVEX(I,2)  =  UG55(I,2) 
UGVEX(I,3)  =  UG55(I,3) 
ELSEIF(N0DE.EQ.  69)  THEN 
UGVEX(I,1)  =  UG69(I,1) 
UGVEX(I,2)  =  UG69(I,2) 
UGVEX(I,3)  =  UG69(I,3) 
END  IF 
605      CONTINUE 
C 

V  =  1 

DO  620  K  =  SMODE, MODAL 

BN(V,1)  =  GAMMA(1,K)*DBLE(UGVEX(K,1)) 
BN(V,2)  =  GAMMA(1,K)*DBLE(UGVEX(K,2)) 
BN(V,3)  =  GAMMA(1,K)*DBLE(UGVEX(K,3)) 
BN(V+1,1)  =  GAMMA(2,K)*DBLE(UGVEX(K,1)) 
BN(V+1,2)  =  GAMMA(2,K)*DBLE(UGVEX(K,2)) 
BN(V+1,3)  =  GAMMA(2,K)*DBLE(UGVEX(K,3)) 
C 

V  =  V+2 
C 

620   CONTINUE 
C 

550   CONTINUE 
C 

C     ***********   ESTABLISH  H,  F  AND  R  MATRICIES    ********** 
C 

DO  50  I  =  1,NH 
DO  55  J  =  1,NH 
H(I,J)  =0.0 
F(I,J)  =  0.  0 
G(I,J)  =0.0 
IF(I.LE.  3)THEN 
L(I,J)  =  0.0 
ENDIF 
55       CONTINUE 
50    CONTINUE 
C 

DO  61  I  =  1,3 
DO  66  J  =  1,3 

R(I,J)  =  0.0 
RINV(I,J)  =0.0 
66       CONTINUE 
61    CONTINUE 
KQ  =  1 
DO  80  K  =  SMODE, EMODE 

H(KQ,KQ)  =  DBLE(LAMA(K)) 
H(KQ+1,KQ+1)  =  1.0D0 
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KQ  =  KQ+2 
80  CONTINUE 

C 

K  =  0 

DO   85   K  =   1,3 

R(K,K)   =  RM 

RINV(K,K)   =   1.  ODO/RM 
85  CONTINUE 

C 

DO   88   I   =   l,2*MODE 

DO  89   J  =   l,2*MODE 
F(I,J)   =  PHI(I,J) 
89  CONTINUE 

88  CONTINUE 

C 

C  *****       COMPUTE  G  MATRIX  AS  NEEDED  BY  RICDSD  SUBR         ***** 

C 

CALL  MATRAN( B , 188 , 2*M0DE , 3 , BT, 3) 
C 

CALL  MATMUL( B , 188 , 2*M0DE , 3 , RINV ,3,3, TEMP ,NH) 
C 

CALL  MATMUL( TEMP , NH , 2*MODE , 3 , BT , 3 , 2*MODE , G , NG ) 
C 

p         y-  y -  j-  y-  .  l.  y-  - ' .  ~  - . .  -  - y  -  y.  y - .  ■  -  y-  -  * . . i  - . .•-  y  -  y-  y - y-  y  -  -  •-  -  •  -  y-  y -  y  -  y- .  •-  -  •  -  y.  y-  -  •  -  y-  y  -  -  •  - y«  y;  -  *  -  y  -  y-  y .  y  -  y  _  -  i -  y.  y  - y  -  y-  y  -  -  *  -  j.  y  -  - '  -  y.  y.  y -  y  -  y. 

C  BEGIN  RICCATI   GAIN  CALCULATIONS 

C»tfi*.M%l*ll*iiSvS»%iS^i»ftrAlfaiV^ 
**   «t    r\   9C   #1    t\   t\   *\    #.    «t   #»    /v    *»    ft    »»    *\    *i    *\    •  -    «i    ^i    '»    *i    '•    #»    ••    r\    *»    *«   4*    ****    •>    •»    '»    *»    '»    •>    /*    '»    '»    (V    *»   *»   #V    '»    r\    *S   *\   t\    »>    »»    '*    '■    /*   -  »    «■    '\    •• 

c 

CALL  RICDSD(NF,NG,NH,NZ,2*M0DE,4*M0DE,F,G,H,Z,W,ER,EI,W0RK, 
+  SCALE, ITYPE,IPVS) 

C 

WRITE(41,*)    '     ' 

WRITE    (41,130) 
130        FORMAT   (/'    THE   CLOSED  LOOP  EIGENVALUES  ARE:  '/) 

DO    140    I   =    1,2*M0DE 

WRITE   (41, *)    ER(I),EI(I) 
140        CONTINUE 

WRITE    (41,150)   WORK(l) 
150        FORMAT  (/'    CONDITION  ESTIMATE    IS: ' ,D26. 18) 
C 

C  ************  COMPUTE   GAIN  MATRIX     L  ************* 

C 

C 

C 


CALL  MATMUL( BT , 3 , 3 , 2*MODE , H , NH , 2*MODE , TEMP 1,3) 
CALL  MATMUL(TEMP1,3,3,2*M0DE,B,188,3,RR,3) 


DO   103   I   =   1,3 

DO    104  J=l,3 

RR(I,J)=R(I,J)+RR(I,J) 
104  CONTINUE 

103  CONTINUE 

C 

CALL  DLINDS(3,RR,3,RRINV,3) 
C 

CALL  MATMUL( RRINV ,3,3, 3 , TEMPI , 3 , 2*MODE , BT, 3) 
C 
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CALL  MATMUL( 61,3,3, 2*M0DE , PHI , 1 88 , 2*M0DE , L , 3 ) 
C 

WRITE(41,*)  '  ' 

WRITE(41,*)  '   GAIN  MATRIX   L  ' 

WRITE(41,*)  '  ' 

WRITE(41,*)  '     ROW  1  ROW  2  ROW  3' 

DO  9155  I=l,2*MODE 

WRITE(41,1040)  (L(J,I),J=1,3) 
9155   CONTINUE 

WRITE(41,*)  '  ' 
C 

Q       *Tfc***Vr****yr*******************Vr*********Vrtfr***:fc**:fc******7V^^ 

C     *****        COMPUTATION  OF  TORQUES  AND  COSTS        ***** 

Q  ********************************************************** 

c 

9000   COUNT  =  0 

TOTCST  =  0.  ODO 

TIME  =  0.  0 
C 

Q  *************************************************************** 

C     *****   SETS  LOOP  FOR  THE  NUMBER  OF  ITERATIONS  NECESSARY   ***** 
C     *****  TO  OBSERVE  THE  SYSTEM  FOR  DESIRED  LENGTH  OF  TIME   ***** 

Q       *************************************************************** 

c 

LOOP  =  INT((MIN*60.  0)/SAMPT) 
PRNT  =  INT(((MIN*60.  0)/SAMPT)/300.  0) 
PRNTG  =  INT(((MIN*60.0)/SAMPT)/1000.0) 
C 

DO  200  N  =0,  LOOP 

TIME  =  DBLE(N)*SAMPT 
C 

IF(N.  EQ.  0)THEN 
IMPLSX  =  IMPX 
IMPLSY  =  IMPY 
IMPLSZ  =  IMPZ 
ELSE 

IMPLSX  =  0.  ODO 
IMPLSY  =  0.  ODO 
IMPLSZ  =  0.  ODO 
END  IF 
C 

Q  ****yrVfVfVffrVfVrVfyf***yfVf***Vf**yfAyryfyf^yfVfyr****************Vf*Vr**VfVfVf** 

C     **********       CONTROL  TORQUE  EQUATIONS       *********** 

C 

SUM1  =  0.  ODO 
SUM2  =  0.  ODO 
SUM3  =  0.  ODO 
C 

DO  210  CT  =  1,  MODE 

CTADJ  =  CT  +  (SMODE  -  1) 

SUM1  =  SUM1  +  L(1,2*CT-1)*X1( CTADJ)  +  L( 1 ,2*CT)*X2( CTADJ) 
SUM2  =  SUM2  +  L(2,2*CT-1)*X1(CTADJ)  +  L(2,2*CT)*X2(CTADJ) 
SUM3  =  SUM3  +  L(3,2*CT-1)*X1(CTADJ)  +  L(3,2*CT)*X2(CTADJ) 
210   CONTINUE 
C 
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c 
c 
c 
c 
c 
c 
c 


c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 
c 


************ 


EVALUATION  AT  ZERO  CONTROL 


************ 


IF(M.  EQ.  1)THEN 

TCX  =  0.  0D0 

TCY  =  0.  0D0 

TCZ  =  0.  0D0 
ELSE 

TCX  =  SUM1*(-1.0D0) 

TCY  =  SUM2*(-1.0D0) 

TCZ  =  SUM3*(-1. 0D0) 
ENDIF 

IF(N.  EQ.  0)THEN 

WRITE  (32,*)  'IMPULSE  X  AXIS,  IMPULSE  Y  AXIS,  IMPULSE  Z  AXIS' 

WRITE  (32,*)  '  ' 

WRITE  (32,1040)  IMPLSX,  IMPLSY,  IMPLSZ 

WRITE  (32,*)  '  ' 

WRITE  (32,*)  'CONTROL  TORQUES  TCX,  TCY,  TCZ' 

WRITE  (32,*)  '  ' 
ENDIF 

IF  (N.  LE.50)  THEN 

WRITE  (32,2000)  TIME,  TCX,  TCY,  TCZ 
ENDIF 

IF  (MOD(N,PRNTG).EQ.  0)  THEN 

WRITE  (32,2000)  TIME,  TCX,  TCY,  TCZ 
ENDIF 


IF  (N.  LE.20 
WRITE ( 
WRITE ( 
ENDIF 

IF  (M.  EQ.  1) 
IF  (MOD(N, 
WRITE(30, 
WRITE(30 
WRITE(31 
ENDIF 
ELSE 
IF  (MOD(N, 
WRITE(30 
WRITE(30 
WRITE(31 
ENDIF 
ENDIF 


)  THEN 

30,1035)  TIME,X1(1),X1(5),X1(10),X1(15),X1(20) 

31,1035)  TIME,X2(1),X2(5),X2(10),X2(15),X2(20) 

THEN 
PRNTG).EQ.  0)  THEN 

1036)  TIME,X1(9),X1(10),X1(11),X1(12),X1(13),X1(14) 
,1036)  TIME,X1(1),X1(4),X1(14),X1(24),X1(34),X1(44) 
,1035)  TIME,X2(1),X2(5),X2(10),X2(15),X2(20) 


PRNT).EQ.  0)  THEN 

,1036)  TIME,X1(1),X1(4),X1(24),X1(44),X1(74),X1(94) 
,1036)  TIME,X1(1),X1(4),X1(14),X1(24),X1(34),X1(44) 
, 1035)  TIME,X2( 1) ,X2(5) ,X2( 10) ,X2( 15) ,X2(20) 


IF  (MOD(N,PRNTG).EQ.  0)  THEN 

COUNT  =  COUNT+1 
ENDIF 

*********************************************************** 

********     SYSTEM  COST  FUNCTION  CALCULATION      ******** 
*********************************************************** 
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SUMC  =  0.  0D0 
ENERGY  =  0.  ODO 
CNTCST  =  0.  ODO 
COST  =  0.  ODO 
C 

DO  230  CF  =  7, MODAL 

MODEN(CF)  =  M0DEN(CF)+(X1(CF)**2)*LAMA(CF)+X2(CF)**2 
SUMC  =  SUMC+(X1(CF)**2)*LAMA(CF)+X2(CF)**2 
230   CONTINUE 
C 

C     ENERGY  =  SUMC 

C     CNTCST  =  (TCX**2)*RM+(TCY**2)*RM+(TCZ**2)*RM 
C     COST  =  ENERGY  +  CNTCST 
C     TOTCST  =  TOTCST  +  COST 
C 

C        IF  (MOD(N,PRNTG).EQ.  0)  THEN 
C  WRITE(33,2000)  TIME, ENERGY, CNTCST, COST 

C        END IF 

IF(N. GE. (LOOP-50))THEN 

WRITE  (34,3002)  MODE, TOTCST 
ENDIF 
C 

Q       *********************************************************** 

C     ********  STATE  UPDATE  EQUATIONS         ********** 

Q       *********************************************************** 

c 

DO  220  KA  =  7, MODAL 
K  =  KA-6 
C 

X1T=PHII(1,1,KA)*X1(KA)+PHII(1,2,KA)^X2(KA)+B((2^K-1),1)^TCX+ 
+       B((2*K-1),2)*TCY+B((2*K-1),3)*TCZ+BN((2*K-1),1)*IMPLSX+ 
+       BN((2*K-1),2)*IMPLSY+BN((2*K-1),3)*IMPLSZ 
C 

X2T=PHII(2,1,KA)^X1(KA)+PHII(2,2,KA)*X2(KA)+B(2^K,1)^TCX+ 
+       B ( 2* K , 2 ) *TC Y+B ( 2*K , 3 ) *TC  Z+BN ( 2*K , 1 ) * I MPLSX+BN ( 2*K , 2 ) * 
+       IMPLSY+BN(2*K,3)*IMPLSZ 
C 

Xl(KA)  =  X1T 
X2(KA)  =  X2T 
C 

220      CONTINUE 
C 

200   CONTINUE 
C 

WRITE  (35,3000)  COUNT 
RTOTAL  =  TOTCST 
C     WRITE  (34,3002)  MODE, RTOTAL 
WRITE  (34,3002)  MODE, TOTCST 
C 

DO  235  K  =  7, MODAL 

RMODEN(K)  =  MODEN(K) 
WRITE  (31,3002)  K,  RMODEN(K) 
235   CONTINUE 
C 

666   CONTINUE 
C 
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* 

* 
* 
* 
* 

* 
* 

* 
* 
* 

* 
* 

* 

* 
* 
* 
* 
* 
* 
* 
* 
* 

J. 

«» 

*505 

C 

C 

C 

C 

C 

700 

701 

702 

703 

704 

705 

706 

707 

708 

1001 

1002 

1004 

1005 

1008 

1010 

1035 

1036 

1040 

1050 

2000 

2001 

3000 


IF(M.EQ.  1)THEN 

MODE=10 

EM0DE=16 
ELSEIF(M.  EQ.  2)THEN 

M0DE=15 

EMODE=21 
ELSEIF(M.  EQ.  2)THEN 
MODE=20 
EMODE=26 
ELSEIF(M.  EQ.  2)THEN 

MODE=25 

EMODE=31 
ELSEIF(M.  EQ.  3)THEN 

MODE=30 

EMODE=36 
ELSEIF(M.  EQ.  4)THEN 

M0DE=35 

EM0DE=41 
IF(M.  EQ.  1)THEN 

M0DE=40 

EM0DE=46 
ELSEIF(M. EQ.  2)THEN 

M0DE=45 

EM0DE=51 
ELSEIF(M.  EQ.  3)THEN 

MODE=50 

EMODE=56 
ENDIF 
CONTINUE 

****************************************  ******************** 

*****  FORMAT  STATEMENTS  ***** 

*********************************************************** 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 

FORMAT( 


( 

( 

( 

( 

( 

( 

( 

( 

( 

IX, 

IX, 

IX, 

IX, 

IX, 

Al) 
t  i 

t  t 

t  t 

i  t 

IX, 

f   t 

14) 


'/STARTING  MODE  NUMBER:   ',12) 
, 'NUMBER  OF  MODES  SCANNED:   ',12) 
, 'NOISE  INPUT  NODE:   ' ,13) 
,' INITIAL  R  VALUE:   *,E12.4) 
,' SAMPLING  TIME:   * ,E12. 4) 
, "DAMPING  FACTOR:   ' ,E12. 4) 
,'LAST  CONTROLLED  MODE:   ',12) 

' /OBSERVATION  TIME:   ',F5.1,'  MINUTES') 

','SIZE  OF  MODAL  MODEL:   ',13,'  MODES') 

A6) 

8E15.  8) 

//) 
60X,E11.5) 

////) 

,F7.2,2X,5(E12.  6,2X)) 

,F7.2,1X,6(E11.5,1X)) 

,3(E15.8,5X)) 

,4(E12.6,2X)) 

F7.2,3X,3(E15.8,3X)) 

,T5,E15.8) 
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3001  FORMAT(F7.2,2X,E12.5) 

3002  F0RMAT(I3,2X,E12.5) 
C 

599   STOP 

END 
C 

C  SUBROUTINE  TO  MATRIX  MULTIPLY 

Q  *y«V}V;toV?Y}Wr?VvV}V}V}V}VycyoViVVfy-VriV:fcy-iV*^^ 

C 

SUBROUTINE  MATMUL(M1 ,LD1 ,R1 ,C1 ,M2,LD2,C2,MP,LD3) 
INTEGER  R1,C1,C2,LD1,LD2,LD3 
REAL*8  M1(LD1,1),M2(LD2,1),MP(LD3,1),SUM 
C 

DO  650  I  =  1,R1 

DO  660  J  =  1,C2 
SUM  =  0. 0D0 
DO  670  K  =  1,C1 

SUM  =  SUM+M1(I,K)*M2(K,J) 
670  CONTINUE 

MP(I,J)  =  SUM 
660  CONTINUE 

650      CONTINUE 
RETURN 
END 
C 

C  SUBROUTINE  TO  TRANSPOSE  A  MATRIX 

C 

SUBROUTINE  MATRAN(MX,LDX,R1,C1 ,MT,LDT) 
INTEGER  R1,C1,IJJ,LDX,LDT 
REAL-'-8  MX(  LDX ,  1 ) ,  MT(  LDT ,  1 ) 
C 

DO  680  I  =  1,R1 

DO  690  J  =  1,C1 

MT(J,I)  =  MX(I,J) 
690  CONTINUE 

680      CONTINUE 
RETURN 
END 
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APPENDIX  B.     KARHUNEN-LOEVE  MODE  PROGRAM 

This  program  computes  the  Karhunen-Loeve  mode  shapes  as  discussed  in  Chapter 
IV.  Subroutine  GENRAA  is  called  to  generate  the  covariance  matrix  of  the  modal  co- 
ordinates, RAA  =  E[a(t)  aT(t)]  .  One  subroutine  from  the  IMSL  math  library 
(DMRRRR)  is  used  for  matrix  multiplication  and  two  EISPACK  subroutines  are  used 
to  compute  the  eigenvalues  and  eigenvectors  of  the  covariance  matrix. 

The  subroutine  GENRAA  as  noted  above  generates  the  covariance  matrix,  RAA, 
of  the  modal  coordinates  of  the  space  station  discussed  in  Chapter  I.   The  algorithm: 

•  Computes  the  natural  frequency,  damping  coefficient,  y,  and  damped  frequency  of 
each  mode 

•  Computes  the  variances  and  covariances  based  on  these  parameters  as  given  by 
equations  19,  20,  21 

The  KL  mode  shapes  are  found  by  solving  the  eigenvalue  problem,  equation  12. 


RAA  *Q  4>  =  /.<j) 


(35) 


where  the  matrix  Q  is: 


Q  =  diag 


C07    0 

|_0       1 

J 

2     t\~ 
(o&     0 

, ...  , 

«56       0 

.  0      l. 

(36) 


and 


•  A  is  the  eigenvalue  of  the  product,  RAA  *  Q 

•  4>  represents  the  orthonormal  eigenvectors 


RAA  *  Q  is  not  a  symmetric  matrix  so  the  eigenvalue  problem  is  changed  to: 

[Q12RAAQ1/2]  [Q120]  =  ;.[Q,20] 


(37) 


Let 


r\  =  Q1  2<£  and  <f>  =  Q"1  2  rj 


(38) 


and  finallv: 


[Q]  2  RAA  Q1  2]  r,  =  ) .n 


(39) 
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The  EISPACK  routines  for  symmetric  matrices  are  used  to  solve  the  eigenvalue  prob- 
lem, equation  39,  and  equation  38  is  used  to  compute  the  KL  mode  shapes. 

This  program  is  run  by  executing  the  following  EXEC  file: 

FORTVS2  KLMODES2 

FI  4  DISK  THESIS  INPUT  B  (PERM 

FI  6  DISK  RAA  OUTPUT  T  (PERM 

FI  8  DISK  EIGVEC  OUTPUT  T  (PERM 

FI  10  DISK  EIGVAL  OUTPUT  T  (PERM 

FI  09  DISK  EIGEX  ERROR  B  (PERM 

FI  12  DISK  XATMODES  D001  B  (PERM 

FI  13  DISK  KLAMAQ  OUTPUT  B  (PERM 

FI  14  DISK  RAA  OUTPUT  (PERM 

FI  17  DISK  NERR001  DATA  A  (PERM 

P  DEF  STOR  2M 

EXEC  MATHPACK  EISPACK 

EXEC  TDISK  5  DIS 

LOAD  KLMODES2  (START 

PROGRAM  KLMODE 
******************  Var  i  ab  1  e  de  f  in  it  ions**  *********************  -k-k-k-k-k-k  *  **** 

*  lama  =  vector  of  natural  frequencies  * 

*  ugvex=  matrix  of  the  modal  amplitudes  and  modal  slopes  * 

*  natmod=vector  containing  the  natural  mode  shapes  that  correspond  * 

*  to  the  KL  mode  shapes  * 

*  klama=vector  of  natural  frequencies  ordered  according  to  th  KL   * 

*  mode  shapes  * 

*  klmod=vector  containing  the  largest  eigenvectors  in  column  of  the  * 

*  eigenvector  matrix  eigenvec  * 

#\  ****/.  I  \     t\     l\     l\     I   \  I   N  M  t\      I   N  1\     *\     .  1  1   *   »  .  /  *  ««  .«  I   *  #*  I*      1\     «f.V  '»  I*     »*  t\     I   \  /.  /  .  1+    t\     I   .  I*     *\  #*  '»  '*     0i    t\    **  "  **  *»  «»  #»  1\    »*  '^    '   *  7*  •   »  "  "  *»  "  #V*t  **  /»  '  *  '»  #.  *»  '  »  *t  #»  *» 

CHARACTER*6  NAM 

DOUBLE  PRECISION  LAMA( 100, 1) ,UGVEX(684, 100) ,D( 100) ,E( 100) 

DOUBLE  PRECISION  KLM0D( 100) , Y,N1 ,Q( 100 , 100) 

DOUBLE  PRECISION  GAMA,DAMP,W(50) ,L(50) ,M(50) ,B2( 1,50) ,Q2( 100, 100) 

DOUBLE  PRECISION  K( 1 ,50) ,A,B,C ,C0EF,RAA( 100 , 100) , APB , APC , TEMPI 

DOUBLE  PRECISION  TEMP2 ,TEMP3,TEMP4,TEMP5 ,TEMP6 ,Q2RAA( 100 , 100) 

DOUBLE  PRECISION  EIGVEC( 100 , 100) ,Z( 100 , 100) ,RAAQ( 100 , 100) 

DOUBLE  PRECISION  Q2RQ2( 100 , 100) ,Q2INV( 100 , 100) ,KLAMA(50) 

INTEGER  KK,LL,ZZ,MODR( 100) ,M0DC( 100) ,MM,NATM0D( 100) ,NUM 

DO  160  1=1,100 

DO  170  J=l,100 

Q2(I,J)=0. 0D00 

Q2INV(I,J)=0.0D00 

RAA(I,J)=0.0D00 
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Q(I,J)=0. ODOO 
170   CONTINUE 
160   CONTINUE 

READ(4,1001)  NAM 

READ(4,1002)  (LAMA( I, 1) ,1=1,100) 

READ(4,1001)  NAM 

READ(4,1002)  ( (UGVEX( I , J) , 1=1,684) , J=l, 100) 
C  NODE  IS  THE  DISTURBANCE  INPUT  LOCATION 
C  AND  DILOC  IS  THE  COMPUTED  LOCATION  OF  THE  MODAL  SLOPE  IN  UGVEX 

NODE=69 

DISLOC  =  4+6*(N0DE-l) 
C  NUMODE  IS  THE  NUMBER  OF  THE  LAST  MODE  CONSIDERED 

NUMODE  =56 

CALL  GENRAA( NUMODE, DISLOC, LAMA, UGVEX, W,L,M,B2,K,RAA) 

WRITE(16,1005)W(1),L(1),M(1),B2(1,1),K(1,1) 

LL=2 

KK=1 

DO  150  3=1, NUMODE 

Q2(KK,KK)=DSQRT(LAMA(J,1)) 

Q2(LL,LL)=1.  ODOO 

Q2INV(KK,KK)=1/Q2(KK,KK) 

Q2INV(LL,LL)  =  1.0D00 

KK=KK+2 

LL=LL+2 
150   CONTINUE 

CALL  DMRRRR( 100 , 100 , Q2 , 100 , 100 , 100 , Q2 , 100 , 100 , 100 ,Q, 100) 

CALL  DMRRRR(100,100,RAA,100,100,100,Q,100,100,100,RAAQ,100) 

DO  15  1=1,100 

WRITE(14,1007)  RAA(I,I) 

WRITE(l.r  ,1007)  Q2(I,I) 
15    CONTINUE 

CALL  DMRRRR( 100 , 100 , Q2 , 100 , 100 , 100 ,RAA, 100 , 100 , 100 , Q2RAA , 100) 

CALL  DMRRRR( 100 , 100 , Q2RAA ,100,100, 100 , Q2 , 100 , 100 , 100 , Q2RQ2 , 100 ) 

CALL  TRED2(100,100,Q2RQ2,D,E,Z) 

CALL  TQL2(100,100,D,E,Z,IERR) 

IF(IERR.NE.  0)THEN 

WRITE(9,*)'ERROR=' ,16 
C     PAUSE 

END  IF 

CALL  DMRRRR(100,100,Q2INV,100,100,100,Z,100,100,100,EIGVEC,100) 

WRITE( 6 , 1005 ) ( (RAA( I , J) ,  J=l , 100) , 1=1 , 100) 

WRITE(8,1006)((EIGVEC(I,J),J=1,100),I=1,100) 

WRITE(10,1007)(D(I),I=1,100) 

MM=1 

DO  10  1=100,1,-1 

KLMOD(MM)=0.0D00 

DO  20  J=l,100 

IF(DABS(KLMOD(MM)).LT. DABS(EIGVEC(J,I)))THEN 

KLMOD(MM)=EIGVEC(J,I) 

MODR(MM)=I 

MODC(MM)=J 

Y=(DBLE(MODC(MM)))/2.  ODOO 

I1=INT(Y) 

N1=Y-(DBLE(I1)) 

IF(DABS(N1).  NE.  0.  0D00)THEN 

NATMOD(MM)=(MODC(MM)+l)/2  +  6 
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ELSE 

NATM0D(MM)=M0DC(MM)/2  +  6 
ENDIF 
ENDIF 
20    CONTINUE 
MM=MM+1 

10  CONTINUE 

DO  30  1=1,100 

WRITE( 12,100)  MODR( I ) ,MODC( I ) ,NATMOD( I ) ,KLMOD( I ) 
30    CONTINUE 

ZZ=50 

DO  180  1=1,100,2 

KLAM A ( Z  Z ) =LAMA ( N ATMOD ( I ) , 1 ) 

ZZ=ZZ-1 
180   CONTINUE 

WRITE( 13 , 1005 ) (KLAMA( I ) , 1=1 , 50) 

NUM=7 

DO  11  1=1,99,2 

IF(DABS(KLMOD(I)).GT.  DABS(KLMOD( 1+1) ) )THEN 

WRITE( 17,200)  NUM,DSQRT(2*(1-DABS(KLM0D(I)))) 

ELSE 

WRITE( 17 ,200)  NUM,DSQRT(2*( 1-DABS(KLM0D( 1+1) ) ) ) 

ENDIF 

NUM=NUM+1 

11  CONTINUE 

100   F0RMAT(1X,I3,2X,I3,2X,I3,2X,E15.  8) 
200   F0RMAT(1X,I2,2X,E15.8) 

1001  F0RMAT(1X,A6) 

1002  F0RMAT(1X,8E15.8) 

1005  F0RMAT(1X,5E15.8//) 

1006  F0RMAT(1X,5E15.  8//) 

1007  F0RMAT(1X,5E15.8//) 
STOP 

END 

SUBROUTINE  GENRAA  (NUMODE,DISLOC ,LAMA,UGVEX,W,L,M,B2,K,RAA) 

DOUBLE  PRECISION  LAMA( 100, 1) ,UGVEX(684, 100) 

DOUBLE  PRECISION  GAMA,DAMP,W(50) ,L(50) ,M(50) ,B2( 1 ,50) 

DOUBLE  PRECISION  K( 1,50) , A, B,C,COEF,RAA( 100, 100) ,APB ,APC, TEMPI 

DOUBLE  PRECISION  TEMP2,TEMP3,TEMP4,TEMP5 ,TEMP6 

DAMP=0. 001D00 

GAMA=DAMP/2.  ODOO 

DO  100  M0DE=7,NUM0DE 

W(M0DE-6)=DSQRT(LAMA  (M0DE,1)) 

L(M0DE-6)=W(M0DE-6)*DSQRT(  1.  ODOO-GAMA*GAMA) 

M(M0DE-6)=GAMA*W(M0DE-6) 
100   CONTINUE 

DO  110  J=7,NUM0DE 
B2(1,J-6)=UGVEX(DISL0C+4*(J-1),J) 
K(l,J-6)=B2(l,J-6)/L(J-6) 
C  K(l,J-6)=UGVEX(DISL0C,J)/L(J-6) 

110        CONTINUE 

DO  130  1=1,50 

DO  140  J=l,50 

A=M(I)+M(J) 

B=L(I)-L(J) 

C=L(I)+L(J) 
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C0EF=K(1,I)*K(1,J) 

APB=2*(A*A+B*B) 

APC=2*(A*A+C*C) 

RAA(2*I-1,2*J-1)=C0EF*(A/APB-A/APC) 

TEMP1=L(J)*B/APB+L(J)*C/APC 

TEMP2=M( J)*A/APB-A*M( J) /APC 

RAA( 2*1-1, 2*J)=C0EF*(TEMP1-TEMP2) 

TEMP3=(L(I)*L(J)*A-L(I)*M(J)*B-L(I)*M(I)*B+M(I)*M(J)*A)/APB 

TEMP4=(L(I)'VL(J)*A-L(I)''-M(J)*C-L(I)*M(I)*C-M(I)'>M(J)*A)/APC 

RAA( 2*1 , 2*J)=C0EF*(TEMP3+TEMP4) 

TEMP5=L(I)*C/APC-L(I)*B/APB 

TEMP6=M(I)*A/APB-M(I)*A/APC 

RAA(2*I,2*J-l)=COEF*(TEMP5-TEMP6) 
140   CONTINUE 
130   CONTINUE 

RETURN 

END 
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